CDOT-CSM-R=-93=-20 


ANALYTICAL SIMULATION 
OF ROCKFALL PREVENTION 
FENCE STRUCTURES 


by 


G.G.W. Mustoe and H.P. Huttelmaier 
Department of Engineering 
Colorado School of Mines 


April 1993 


Prepared under contract with the 
Colorado Department of Transportation 
in cooperation with the U.S. Department 
of Transportation, Federal Highway 
Administration 


Technical Report Documentation Page 


1. Report No. 2. Government Accession No. 3. Recipient’s Catalog No. 
CDOT-CSM-R-93-20 


4. Title and Subtitle 5. Report Date 
Analytical Simulation of Rockfall Prevention Fence April 1993 
Structures 6. Performing Organization Code 


G.G.W. Mustoe and H.P. Huttelmaier CDOT-CSM-R-93-20 

9. Performing Organization Name and Address 
Colorado School of Mines 
Department of Engineering 
Golden, Colorado 80401 

12. Sponsoring Agency Name and Address 
Colorado Department of Transportation Epo 
4201 East Arkansas Avenue 14. Sponsoring Agency Code 
Denver, Colorado 80222 oa 


15. Supplementary Notes 


Prepared in Cooperation with the U.S. Department of Transportation Federal 
Highway Administration 


16. Abstract 

This report describes the research work performed for the Colorado Department of 
Transportation (CDOT) during the period from June 1991 to March 1993, under the CSM 
Contract No. 3422. 

The work performed during this contract included the following tasks: (i) the 
development of a novel numerical model for the dynamic simulation of a rockfall prevention 
fence using the discrete element method (DEM), (ii) the numerical implementation of the 
DEM numerical model in a computer code to be used by engineers at CDOT, and (111) 
demonstration and validation of the DEM computer technique to simulate rockfall impacts 
against rockfall fences. 

The original DEM numerical mode] was designed to analyze flexible rockfall prevention 
fences that were comprised of a single layer of columnar attenuators. In an additional 
research task funded under this contract the DEM model was extended to deal with fences 
with two layers of attenuators. 


17. Key Words 18. Distribution Statement 
Rockfall Prevention Rockfall Structures No Restrictions: This report is 


Discrete Element available to the public through 
Method | the National Technical Info. 
Service. Springfield, VA 22161 


19.Security Classif. (report) 20.Security Classif. (page) 21. No. of Pages 
Unclassified Unclassified 79 


1 


TABLE OF CONTENTS 


Page No 
SUMMARY. jscscsscctcte tee Wsscade cel tecaeitie et hicaateceteel er uladen te Moped necslesaun ita 4 
LO INTRODUCTION 2h ccorccredichscheciuauaei active A ea deian eames a er aaeamusendnmeaeaions 4 
1.1 The Prototype Fence Structure and Impacting Rocks ..................: > 
1.2 The Discrete Element Method ...............:cccecsessseceeseeeseeeeeeesceeeescaeens 5 
20) PAE PENCE MODEL coset chesedcssticeetenienaveiee wecidnn sere teeet cee eeedese 6 
2.1 Rigid Body Dynamics of the Attenuator Columns and Rocks ........... 6 
2.2 Rigid Body Geometry y..c.c:sisschasareecelarinte inden Adeceae 6 
2.3 Impact Geometry Conditions) ....:.sescecccccsiestencontsteaterececsnsersscsdenceeatoustes y | 
2A Cable Connect Ons apc .tidadetcsaragiar tan dete wlevcteectadatdtteactamsaciotereheanes 7 
2.5 Dynamic Contact Forces sa... 1.<./scc¢ececkcstates teri a lehaorasteceaataebedaanaeeentns 7 
3.0 COMPUTATIONAL DETAILS OF FENCE MODEL ...000.....eseeceeeseeeeeeteeees 8 
3.1 Attenuator Equations of Motion ...............:cc:ccsceceseeseeeeresrceceteesreneeeees 8 
3.2 Rock Equations Of Motion ................ccsccccssesescenecsceseseeceetaeceecteetseeeaes 9 
3.3 Top Cable Attachment Forces ..0........ecec eects esaseeeesesesaseaeseesenseesens 9 
3.4 Bottom Cable Attachment Forces ...............::ccccscscceseceseeseeeseetesaeeeers 10 
BD COMACE POLES cps cencsct cpacepenasvexeacataioneaecucanariis eieiesaveme ee ainkoeeas 11 
3.6: Restitution Model 5.4.4. ce vcesccr cue epncds eae dense ovieetbaaeians 12 
3.7 Time Steppins Algorithm i331 ds ee iiscsccdeeooceieedaeceessidaveees 13 
4.0 DESCRIPTION OF THE DEM SOFTWARE ..............::ccscsscsssssssceeesesseneeees 14 
5.0. VALIDATION. PROBEEMSG .auz.cessegcsscaenavapaosncagsnpeseasassonetaoeeooueatesonvoratiaies 15 
60; ROCK IMPACT STUDIES =. cesatcGuetctaa a aie ene 17 
6.1 Model Dimensions and Properties ..............:.c::csccesccsnessseseresnsesceeers 17 
6.2 Sunilation, Results: gs. cocesaceususeaaresngss aaewsare tected meee semen 18 
O:3 ODSERVALONS cpesccescrecteas Antes seid aoe alent ete nen tee 18 
7,0 CONCLUDING REMARKS sies.cicceicereeceustanatettece stead cs aacceieeesentadm cies 19 
SO AKNOWLEDGEMENTS scisdsicsechatesielalsuganvlasnnlel decent iaaian luis. 20 
OO REFERENCES 55 cSccs, cease lac ctendsenceteate cassette ctr aneetn Maier ee 20 
APPENDIX I - Modeling of Single and Double Layered Flexible Rockfall 
Prevention Fences <5. ca.uiriy divas Creates ths ral ee ede hace one te taeita El 


APPENDIX II - Computer Instructions and Users Manual for FENCE 


SUMMARY 


This report describes the research work performed for the Colorado Department of 
Highways during the period from June 1991 to March 1993, under the research contract 
entitled Analytical Simulation of Rockfall Prevention Fence Structures, CSM Contract 
No. 3422. 


The work performed during this contract included the following tasks:(i) the development 
of a novel numerical model for the dynamic simulation of a rockfall prevention fence using 
the discrete element method (DEM), (ii) the numerical implementation of the DEM 
numerical model in a computer code to be used by engineers at the Colorado Department 
of Highways, and (iii) demonstration and validation of the DEM computer technique to 
simulate rockfall impacts against rockfall fences. 


The original DEM numerical model was designed to analyze flexible rockfall prevention 
fences that were comprised of a single layer of columnar attenautors. In an additional 
research task funded under this contract the DEM model was extended to deal with fences 
with two layers of attenuators. The details of this additional research work are contained 
in Appendix I. 


The dynamics simulation software was developed on an Intel based CPU 80486 
microcomputer. This hardware environment was also used to perform the dynamic 
analyses which are described herein. 


1.0 INTRODUCTION 


The purpose of a rockfall prevention fence is to decrease the kinetic energy of a falling 
rock or rocks down a mountain valley. Such fences are installed near mountain highways 
to minimize the risk of automobile accidents caused by rockfall events. The fence type 
being investigated consists of hanging columnar attenuator masses assembled of used 
truck tires and rims, threaded onto a steel rod, which are suspended from a horizontal 
overhead wire cable and connected by a bottom cable (see Fig. 1). When a rock hits the 
fence the hanging columns begin to move in a pendulum-like motion. This mechanism 
transfers a significant amount of the rock's kinetic energy during the impact phase. Since 
the fence structure is flexible and not rigid, the damage to the fence during most rock 
impacts is only slight. In many cases the decrease in rock velocity after its impact with the 
fence is sufficient so that the rock will harmlessly come to rest shortly after impacting the 
fence, and not enter the highway. In extreme cases where very large rocks (over 2000 Ibs) 
moving with high velocity hit the fence, the fence capacity is exceeded. In this situation the 
hanger connection at the top of the steel rod that passes through the center of attenuator is 
designed to break. This failure results in the loss of an attenuator column but not the 


complete breakdown of the fence structure. This design feature allows the fence to be 
repaired easily, with the simple reattachment of a new attenuator column [2]. 


The numerical model reported herein, provides a technique that can be used in 
conjunction with experimental methods to estimate the fence capacity. The present study 
determines the effects of varying fence design and rock impact parameters on the 
magnitudes of the relevant dynamics response parameters. For example, the dynamic 
impact calculations study the effects of different rock sizes and rock velocities on 
particular fence configurations. Fence response parameters of particular importance 
include; (i) the forces in the hanger connection and the overhead and bottom cables, and 
(ii) the decrease in rock velocity due to impact. 


The numerical procedure developed to perform this dynamic impact simulation is based 
upon the DEM. This procedure solves the nonlinear coupled dynamic equations of motion 
for the idealized fence structure and the impacting rock. The DEM model assumes that the 
fence structure can be idealized as a system of connected rigid bodies, and the rock as an 
impacting rigid body. The dynamic fence-rock impact conditions are modeled with a 
penalty formulation and an automatic contact detection algorithm. This numerical 
procedure utilizes an explicit time stepping scheme to solve the discretized equations of 
motion. 


1.1 The Prototype Fence Structure and Impacting Rocks 


A typical fence structure, similar to those presently being investigated, is shown in Fig. 2. 
Cylindrical column attenuator masses are connected to a steel wire overhead cable through 
hanger attachment cables. These column masses consist of layers of used truck tires on 
their rims that are threaded onto a vertical steel rod. The introduction of the bottom cable 
connection ensures that the combined mass of all attenuators is active during an impact. It 
should be noted that any slack in the bottom cable results in a delayed displacement action 
between attenuators. The area of the bottom cable is approximately one quarter of the 
area of the main overhead cable with the same modulus of elasticity as the overhead cable. 
The bases of the attenuator columns are close to the ground to capture low falling rocks. 
In some cases the attenuator columns are of different lengths to adapt to the uneven 
profile of the terrain. 


The sizes of the most likely impacting rocks range between 500 Ibs to 1500 lbs. 
Corresponding impact rock velocities vary between 50 and 70 ft/s. It should be noted that 
larger impacting rocks that are over 2000 Ibs in weight are likely to cause significant fence 
damage, and exceed the fence capacity. 


1.2 The Discrete Element Method 
Discrete element methods are a family of related numerical techniques specifically 


designed to solve problems in applied mechanics, which exhibit gross discontinuous, 
material and geometrical, behavior. For example, many DEM's are used to analyze systems 


of interacting rigid or deformable bodies undergoing large dynamic or pseudo static 
motion, governed by complex constitutive behavior. It should be noted that the solutions 
to problems of this type are often intractable by conventional continuum based methods 
such as finite element, finite difference and boundary element procedures. 


A typical discrete element algorithm includes the following features: 


Gi) an automatic contact detection scheme which can recognize new contacts as the 
simulation progresses, and update the contact topology appropriately, 

(ii) a contact force model to determine the forces acting between interacting bodies, and 
(iii) a time integration scheme to solve the governing equations of motion, which allows 
for finite displacement and rotation of the discrete bodies within the system. 


Applications of DEM's include; the mechanical modeling of granular media, failure 
analysis of brittle materials (ceramics, rocks, ice, etc.), mechanical behavior of fractured 
and jointed rock masses, and discrete simulation of compaction in material processing of 
ceramics and powder technologies [3,6]. More recently, discrete elements have also been 
applied to continuum structural mechanics problems, which include large displacement 
dynamic analysis of elasto-plastic beams [4], and deep-ocean mining pipes [5]. 


2.0 THE FENCE MODEL 


The modeling aspects and assumptions of the discrete element model, namely the rigid 
body dynamics of the fence and rock, the cable behavior, the rock/attenuator impact 
interaction model and the numerical integration of the dynamic equations of motion are 
discussed below. 


2.1 Rigid Body Dynamics of Attenuator Columns and Rock 


The attenuators and the impacting rock are both represented by three-dimensional rigid 
bodies. The local elastic or plastic deformation that occurs in the vicinity of the point of 
contact during the impact is accounted for within the contact force law. Details of the 
contact model are given in a subsequent section of this paper. The fence attenuator 
columns are modeled as circular cylinders, and the impacting rock as a sphere. The 
kinematics of the fence attenuator columns and the rock is assumed to be planar. The 
dynamic motion of the fence attenuators and rock are therefore, confined to the x, y plane 
(see Fig. 2b). 


2.2 Rigid Body Boundary Geometry 
A superquadric geometrical formulation is used to describe the two- dimensional surface 


geometry of the outline of attenuators and the rock. This geometrical surface model can 
be used to describe any shape, from an ellipse to a rectangle. The boundary shape of the 


two-dimensional body is defined by a number of boundary nodal points on the super 
ellipse: 


(x/a)” + (y/b)" =1 (1) 


Notes: (i) x, y are principal local centroidal coordinates for the super ellipse, and (ii) the 
parameters a, b and n define the bounding dimensions and shape of the super ellipse. 


The current positions of these boundary nodes are then utilized to check for geometrical 
contact or overlap with neighboring bodies. This check is performed with a simple 
algebraic calculation of an inside-outside function. It should be noted that the superquadric 
representation provides analytical formulae for the convenient computation of the outward 
normal and curvature at any point on the boundary of the body. For further details, the 
reader is referred to [1]. 


2.3 Impact Geometry Conditions 


It is assumed that the rock hits the central portion of the fence (Fig. 1). This assumption 
allows the use of symmetry. Furthermore, since the average dimension of the rock is 
approximately equal to the diameter of an attenuator, if the fence contains an odd number 
of attenuators, the rock impacts the central attenuator only, otherwise, the rock impacts 
the two central attenuators symmetrically and the impact loads are equally distributed 
between the two central attenuators. In this situation the two central attenuators are 
assumed to move together. 


2.4Cable Connections 


The sections of cable connecting adjacent attenuators are idealized as massless cable 
elements with linear elastic or perfect elastic-plastic material behavior. The pre-impact 
position of the overhead cable is defined as the static equilibrium shape due to the self 
weight of the attenuator columns. This position is simply computed by a direct application 
of the static equilibrium equations for the fence. These cable elements transmit tension 
only and may become slack if appropriate kinematic conditions occur. Note, that since the 
attenuators are constrained to move only in the x, y (Fig. 2b) plane, the z-coordinates of 
the ends of the cable elements do not change. 


2.5 Dynamic Contact Forces 


The dynamic contact forces generated during the impact are determined with a penalty 
formulation. The penalty method employs a contact "element" at the point of impact 
between the two bodies that is defined by a combined concentrated stiffness and damping 
element. This contact "element" produces forces that depend upon the current penetration 
or overlap distance and relative impact velocity between the impacting rock and the 
attenuator at the point of impact. These forces are applied in an equal and opposite 


manner to the rock and the attenuator. Further details of the penalty approach are given in 
a subsequent section of this report. 


Note, that only the central attenuator(s) is impacted. The other attenuators act as slaves 
that move only because of forces transmitted through cables connecting adjacent 
attenuators. Furthermore, note that contact forces between adjacent attenuator columns 
are neglected, since they are small compared with the impact force between the central 
attenuator and the impacting rock. 


3.0 COMPUTATIONAL DETAILS OF FENCE MODEL 


3.1 Attenuator Equations of Motion 


The equations of motion for the ith fence attenuator, assuming that the dynamic motion is 
planar and confined in the x-y plane are given by: 


m,; xX; = 2, (2) 
m; 9; = ZF, (3) 
1,6, ==M, (4) 


for i= 1,2,...N, where N, is the total number of fence attenuators, and N is defined as N 
= (N, + 1/2. See Fig. 3 for details of the attenuator identification numbering scheme. The 
corresponding resultant loading terms in equations (2),(3) and (4), are: 


LF, = 1, +B, +P, (5) 


LF, =1, +B,+P,-mg (6) 


=M; = B,,bcos0; +B, bsin®, + T,acos@, + T, asin®, 
+P, Ay; + P, Ax; (7) 


where T., and T, are the global x- and y- components of tension forces induced at the top 
overhead cable connection, B,, and B, are the global x- and y- components of tension 


forces induced at the bottom cable connection, P,, and P, are the global x- and y- 
components of the contact forces generated by the rock attenuator impact, and g is the 


gravitational constant. For the ith attenuator the following parameters are defined: m, is 
the mass, I, is the mass polar moment of inertia with respect to the centroid, g is the 
gravitational constant, x, and y, are the global centroidal coordinates, and CJ is an angular 
coordinate defining the axial direction of the attenuator (see Fig. 4). The distances of the 


top and bottom cable connections with respect to the attenuator centroid are defined by a 
and b, respectively. The global x- and y- coordinates of the point of contact of the rock 


with the attenuator measured with respect to the attenuator centroid are defined by Ax, 
and Ay,. Note, that if N is even, the mass and moment of inertia terms in equations (2) 


through (4), for i=1, are doubled. This modification is required because of the assumed 
symmetrical impact geometry. 
3.2 Rock Equations of Motion 


The corresponding equations of motion for the rock are: 


m,x, =P, (8) 
m, Y; = Py —m, g (9) 
1,6, =P, Ax, -P, Ay, (10) 


where m, is the mass of the rock, I, is the rock mass moment of inertia with respect to its 


centroid, and Ax, ,Ay, are the global x- and y- coordinates of the point of contact of the 
rock with the attenuator, with respect to the centroid of the rock. 


3.3 Top Cable Attachment Forces 


In Fig. 5, the top overhead cable segment with end points T, andT,,, between attenuators 
1 and i+] is shown. In a coordinate system with the origin at point T,, the coordinates of 


pointT,,, are denoted by (x,,,, Yis,, d;), where d, is the horizontal distance between 
attenuators i and i+1, which is constant because of the assumption of planar kinematics. At 
time t = 0, the initial position of the pointT,,, with respect to point T, is given by: 


(9, Yin, —din) = (0, y;, 0) + (0, H, -di,,) (11) 


where H=h,,, — h,. Note that the initial vertical positions,h,, which define the pre- 
impact configuration are calculated from static equilibrium as discussed previously. The 
initial unstretched length of the cable segment at time t = 0 is computed from: 


£4 = (dj +H’) -T,/k, (12) 
where k, is the axial stiffness of the cable segment, and T,, is the initial tension of the 
cable segment in the initial equilibrium condition. The current length of the cable segment 
at time t is given by: 

€,=(x¢ ty; tdi)” (13) 


The tension in the cable segment at time t is computed from: 


(14) 


i 


aon k(@;-£,) for €; 22, 
Bie’ for £; <2, 


Note, that the above equation models cable slack correctly and non-physical compressive 


forces cannot be transmitted. The corresponding global x- any y- components of the cable 
tension are given respectively by: 


T, =T,x,/@, and Ty, =Tyy,/2; (15) 


The resultants global x- and y- components of forces acting on the ith attenuator due to 
the tension at the top cable connection induced in the cable segments between the (i-1)th, 
ith and (i+ 1)th attenuators is determined from: 


%=]%—-T4. and T, =T, -T (16) 
It should be noted that for the first attenuator (i=1),T,, =T,, =0, and because of 


symmetry the resultant forces T,, and T), are doubled. 


3.4 Bottom Cable Attachment Forces 


The previous analysis for the computation of the top overhead cable forces is applicable to 
the bottom cable forces with the following minor modifications: 


(1) the axial cable stiffness for the bottom cables should be changed accordingly, and 
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(ii) the initial tension in the bottom cables is zero, and there is a significant amount of slack 
cable length which must be taken up before any tensile forces are transmitted to the 
attenuators. 


Condition (i) requires that the axial stiffness of the bottom cable is doubled because of the 
double cable wrapping in the fence design (see Fig. 2a). 


Condition (ii) is accounted for by modifying the corresponding form of the equation (12) 


for the unstretched cable length by setting T,, = 0 and by increasing the unstretched length 
by a user-defined slack cable length. 


3.5 Contact Forces 


An accurate representation of the dynamic impact forces, stresses and strains generated 
during the fence-rock impact event would require a detailed non-linear dynamic impact 
analysis. This could be performed with a finite element calculation in which the rock and 
the attenuators would be idealized as deformable bodies. However, since only the energy 
transfer mechanism between the rock and the fence is of primary interest, a simple 
dynamic contact model can be used. The normal contact element model is comprised of a 
linear spring stiffness and viscous damper, connected in parallel, that acts in the contact 
normal direction at the instantaneous point of contact, as shown in Fig. 6a. Behavior in the 
contact shear direction is defined by a shear contact element model consisting of a linear 
spring stiffness and a slider that simulates a simple Coulomb dynamic friction model (see 
Fig. 6b). Note, that the contact normal and shear directions are defined by the boundary 
geometry of the body with the smoothest tangent plane at the point of contact. The normal 
and shear contact element model are active only when the discretized boundaries of the 
rock and the attenuator penetrate each other. 


The normal instantaneous contact force, F,, is defined by: 


EK, = (17) 


n 


k, 8, +¢8, for 5,20 
0 for 5, <0 


where, k, is the normal contact stiffness, c is the normal contact viscous damping 
coefficient, and 5, is the normal penetration distance at the point of contact. 


The shear instantaneous contact force, F., is given by: 


k,5, for k,|8,|<w|F| and 8, >0 
F={pF, for k,|6,|2u|F| and8, >0 (18) 
0 for 5, <0 


where, k, is the shear contact stiffness, 5, is the shear slip distance at the point of contact 
and pW is the dynamic coefficient of friction between the rock and the attenuator. 


Note: A suitable value of the normal contact stiffness k, is specified by the user by 
limiting the maximum penetration, or overlap, between the surfaces of the impacting two 
bodies to a predefined value. Typically, this maximum penetration is defined as a small 
fraction of the average body dimension of the contacting bodies. For example, this fraction 
is usually between 0.01 and 0.05. An estimate of k, can simply be determined with an 
approximate energy calculation for the impacting bodies. 


For fence-rock impact studies the normal contact stiffness k, should be determined by 
prescribing an approximate maximum overlap between the rock and attenuator elements. 
A facility for this has been included in the DEM computer model and is mentioned in 
Appendix II of this report. In the analyses performed in Appendix I for double layer fences 
an overlap of six inches was prescribed. 


3.6 Restitution Model 


In the discrete element numerical model, the energy loss due to impact ts simulated with a 
viscous damper which acts at the point of impact, P, in the contact normal direction. The 


normal viscous force, F,, in this damper, which is the second term in the normal contact 
force equation (17), is given by: 


F=cd (19) 


where 5, is the normal component of the relative velocity of the fence with respect to the 
rock at the point of contact, and is defined by: 


5, = (Voe+ Vo) i (20) 


Note, that 0 is the unit outward contact normal direction at the point of impact, P, and V,, 


and V,, are the velocities of the fence and the rock during the impact, respectively, and c is 
the normal contact viscous damping coefficient as defined in equation (17). 


In a dynamic impact between two bodies, the energy loss can be approximately modeled 
with a simple restitution model which requires e, the coefficient of restitution, to be known 
a priori via experimental data. Due to lack of available quantitative data concerning the 
impact energy losses during a fence-rock impact situation this simple model is appropriate. 
However, in order to apply the normal contact force law, equation (17), a relationship 
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between the viscous damping coefficient c, and the coefficient of restitution e, is required. 
For this purpose, a simple approximate relationship between the viscous damping 
coefficient, c, and the coefficient of restitution, e, can be derived by assuming that the rock 


and the attenuator are approximated as two particles of mass m, and m,, respectively. 
This results in the relationship: 


c= 2In(1/e)[k, m,m,/(m,+m,)]"’/ [n? +(In(1/e))?]” (21) 


where k, is the normal contact stiffhess at P, the point of impact. 


It should be noted that the use of equation (21) is only approximate when applied to a 
non-central impact between a spherical rock and a cylindrical attenuator. Some numerical 
experimentation, therefore, is usually required to obtain the precise value of the viscous 
damping coefficient which corresponds to a specified value of the coefficient of restitution. 


3.7 Time Stepping Algorithm 


The above equations of motion, (2), (3) and (4), which describe; (i) the translational rigid 
body motion of the center of mass for each attenuator element, and, (ii) the rotational 
motion of each element, can be written in the generalized form: 


My i =EX, k=1,2,3 (22) 


where q; is a generalized coordinate, My is a generalized mass, and F* is a generalized 
loading vector. 


For a planar rigid body element, the generalized variables are identified as follows: 

(i) the generalized coordinates, q;, q?, are the element centroid Cartesian coordinates x,, 
y,, and q; is the angular coordinate @,, of the elements longitudinal or local x-axis, 

(ii) the generalized masses, M;, M; are the element mass m,, 

and M? is the elements polar mass moment of inertia I,, and 


(iii) the generalized loading vector, F’, F’, are the x,y- components of applied forces 
acting at the element centroid and 


F° is the applied moment acting at the element centroid. 
The solution of the above decoupled dynamic element equations for large rigid body 


motions is obtained with an updated Lagrangian algorithm. This procedure employs 
incremental updating of the element centroid kinematic variables. Time integration of the 
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dynamic element equations is then performed with an explicit central difference time 
stepping scheme. 


The_central difference time step update scheme for the generic dynamic equilibrium 
equation (22), from time t* to t” + At, is performed in two parts: 

(i) a generalized velocity update for the linear and angular velocity of the center of mass of 
all the elements, where the generalized velocities, qx, for a typical element are updated by: 


mti/2 qk = n2gk + At® FE /Mé (23) 


1 


where the pre-superscripts n-1/2, n, and n+1/2 indicate quantities defined at times 
t” —At/2, t", and t" + At/2, respectively, and 


(ii) a generalized position update for the coordinates of the center of mass and the local 
longitudinal axis of all the elements, where the generalized coordinate for a typical 


element, qi, is updated by: 
mil gk = "qi + mza* At (24) 


The pre-superscripts n and n+] indicate quantities defined at times t" and t” + At. 


The numerical stability criteria of this time stepping scheme depends on the maximum 
frequency of the combined discrete element system, w,,.. This gives a time step limitation 
constraint of the form: 


At nax < 2! Drre (25) 


where At,,,, is the maximum allowable time step. This criteria can be applied simply with 
the largest natural frequency of any rigid attenuator element within the fence structure 
which is an approximate upper bound of @,,,. This approximate procedure avoids a 
costly eigenvalue analysis of the whole discrete element fence idealization. 


4.0 DESCRIPTION OF THE DEM SOFTWARE MODULES 


The DEM dynamic fence-rock analysis software consists of three major modules: 
(i) the computational DEM module (FENCE), 

(ii)the geometrical graphics snapshot and animation module (GPLOT), and 
(iii)the time history graphics module (HPLOT). 
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In the DEM computational module, FENCE, the input data contains the following 
parameters: 


@) fence and rock geometry and mass properties, 

(ii) material model and properties for the fence cables, 

(iii) impact and restitution parameters between the fence and the rock, 

(iv) the initial rock velocity prior to impact, and 

(v) control parameters which describe the type and time intervals of the computed results. 


This input data is contained in a data file with a user-specified name appended with the 
extension .DAT. 


The computed results are written to three different output data files with a user-specified 
name appended with the extensions .OUT, .HST, and .GEO, respectively. These output 
data files are called the printout file, the time history file, and the geometry file, 
respectively. 


The printout file is used directly to verify the input data and for checking purposes by the 
user. The time history and the geometry results data files are input data for the graphical 
post-processing modules HPLOT and GPLOT, respectively. 


The HPLOT graphics module interactively displays time history plots and produces 
postscript graphics output. Typical time history data plots include; (i) centroidal 
kinematics quantities such as position and velocity components for the rock and the 
attenuator columns, and (ii) forces in the overhead and bottom cable segments, and the 
attenuator hanger connections. 


The GPLOT graphics module is used to display the following geometrical data for the 
fence-rock system; (i) single snapshots of the fence-rock geometry at prescribed times 
during the impact, or (ii) animation of the fence-rock motion throughout the impact. 
Postscript graphics output of the single geometry snapshots can also be produced by 
GPLOT. 


The hard copy graphics output from HPLOT and GPLOT have been used to illustrate the 
computational results obtained from the numerical analyses presented in the subsequent 
sections in this paper. It is the opinion of the authors that without graphical output these 
analyses would be very difficult to perform and interpret effectively. 


5.0 VALIDATION PROBLEMS 
The following example problems are presented to validate the above described discrete 


element computational procedure. These calculations illustrate that this numerical 
approach accurately models; (i) dynamic impact phenomena between rigid bodies, (ii) 


large dynamic rigid body motion and translation, and (iii) large elastic deformations of the 
flexible cables attached to the bodies. 


The first example problem consists of a particle of mass m, initially at rest, attached to two 
unstretched elastic cables of stiffness, k, and length, 2, in a horizontal plane (note that the 
gravitational effects are ignored), being impacted by a particle of mass M with an initial 
velocity Vo (see Fig. 7). The coefficient of restitution for this impact is defined by e. 


The data used for this analysis are as follows: 
k=4N/n, € =10m, m=1kg,M=2kg, e=0.5, and vo =2 m/s. 


The analytical solution for this problem is described by: 
v, = M(i+e)v, /(M+m) (26) 
Where v, is the velocity of the particle of mass m, immediately after the impact, and 


v?(8)= v2 —2k é?(1—-cos@)’ /(m cos? 9) (27) 


where v(@) is the velocity of the particle of mass m after the impact at the instant when 
the cables subtend an angle @ with respect to their initial directions. 


The time histories of the x-component of velocity for the particles of masses m and M are 
shown in Fig. 8. The distance moved by the particle of mass m, after its velocity decreases 
to zero is computed as 3.826 m and corresponds to an angle @ = 0.365 radians. This is in 
good agreement with the exact result which can be obtained ffom equation (27). It should 
also be noted that the velocity immediately after impact is 1.997 m/s, which compares 
closely with the analytical value of 2.0 m/s, predicted by equation (26). This result 
validates the viscous damper - restitution coefficient equation stated in equation (21), 
previously. 


The second validation problem consists of a single cylindrical attenuator being impacted 
by a spherical rock moving with an initial horizontal velocity vo. The cylindrical 
attenuator is pin-jointed at its top, and the line of impact is through the center of mass of 
the attenuator, (see Fig. 9). The superquadric geometrical representation of two- 
dimensional outline shapes of the spherical rock and the cylindrical attenuators in the 
subsequent analyses are defined as follows: 

(i) the rock geometry is defined by a circular boundary shape with twelve equally spaced 
nodes, and 

(ii) the fence attenuator geometry is defined by a boundary shape with four nodes at the 
comers of the rectangular boundary. 
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The attenuators mass, moment of inertia with respect to the center of mass, vertical height 
and radius are defined by m, Ig, @ and r, respectively. The mass of the rock is M, and the 
coefficient of restitution is e. The data for this analysis are given as: 


m = 14.92 slugs, 2 = 8 ft, Ig = 85.42 slug-ft2, r= 1.2 ft, 
M = 37.23 slugs, Vg = 5 ft/s and e= 0.5. 


The analytical solution for this problem is given by: 


v,=M(1+e)(¢/4)/[I,+(M+m)é?/4] (28) 


where v1 is the speed of the center of mass of the attenuator immediately after the impact. 
Note: the corresponding angular velocity is determined by: 


@,=2v,/£ (29) 


The subsequent angular velocity of the center of mass of the attenuator, @(@), when the 
angle of the longitudinal axis of the attenuator with respect to the vertical direction is 0, is 
defined by: 


(0) = w? + 2mgé(cosO-1)/(I, +m 2? / 4) (30) 


The numerical results of the dynamic impact simulation are given in Figs. 10 and 11. Fig. 
10 shows snapshots of the geometry of the rock and the attenuator during the dynamic 
impact simulation. A time history of the horizontal component of velocity of the rock and 
the attenuator during and after the impact is shown in Fig. 11. The computed horizontal 
components of the rock and the attenuator immediately after the impact are 2.36 ft/s and 
4.85 fi/s respectively, which are in close agreement with the exact solution. The maximum 
angle of rotation of the attenuator after the impact was computed to be 0.503 radians 
which agrees with the analytical value to three significant figures. This analytical result can 
be computed from equation (30). 


6.0 ROCK IMPACT STUDIES 


6.1 Model Dimensions and Properties 


A full size fence model is used to study the dynamic response of the fence due to 
impacting rocks. A similar full size prototype fence has been built by the Colorado 
Department of Highways (CDOH). Limited experimental data, which is mainly qualitative, 
has been obtained from this fence structure by the CDOH. Dimensions and material 
properties for this numerical idealization are described in Fig. 2 and Table 1. In the present 
computer model two idealized material models for the cable behavior have been 


implemented. The first is a simple linear elastic model, and the second is a perfect elastic- 
plastic model. These models were selected, since more precise cable data for the 
prototype fence was not available. It should be noted that more realistic material models 
for the cable can be readily inserted into the present algorithm. For example, a 
constititutive model with time dependent stiffness properties which accounts for cable life 
[7], could be implemented. 


6.2 Simulation Results 


A series of analyses were performed which included the following possible combinations: 
(i) two rock sizes, (ii) the two material models, (iii) different rock impact heights, (iv) 
different bottom cable slack conditions, and (v) no bottom cable. Each computer run 
required approximately 100,000 time steps with At = 2 x 10-6 seconds. Each analysis 
required about 15 minutes, using an Intel CPU based 80486, 33 MHz microcomputer. 


Tables 2 and 3 summarize the computed results for the 600 lb. rock, and the 1200 lb. 
rock, respectively. The calculations include the maximum forces in the fence cable 
segments and the attenuator hanger connections, and the final rock velocities after the 
impact. The results quoted are for analyses performed with; (i) the linear elastic cable 
model, and (ii) the perfect elastic-plastic cable model. Note, that the results for the perfect 
elastic-plastic model are quoted in parenthesis in Tables 2 and 3. A bottom cable slack of 
0.5 ft. is assumed. The first column in Table 2 indicates the vertical variation of the rock 
impact position with respect to the attenuator centroid. 


Figures 12 through 15 show detailed geometry snapshots and time history plots for an 
impact situation, when the rock impact location coincides with the centroid of the 
attenuator. Figs. 12 and 13 show snapshots of the fence-rock geometry during the impact 
for the 600 Ib rock, with a bottom cable (0.5 ft slack), and without bottom cable, 
respectively. Figures 14 and 15 illustrate time histories for the primary design parameters, 
namely, attachment force, forces in top and bottom cable, and the rock horizontal velocity. 
These results were computed for both rock sizes and cable material models. The 
computed results presented in Tables 2 and 3, and the time history data in Figs. 14 and 15 
are discussed below. 


A series of further analyses were performed for 600 Ib and 1200 Ib rock impacts against a 
single layer fence structure. In these simulations the bottom fence cable was analyzed 
under the two following conditions: (i) 0.0 ft. slack, and (ii) 0.25 ft. slack The maximum 
forces in the top overhead cable and top joint at the central attenuator are quoted in 
Tables 4 and 5. From the results of the analyses it was concluded that different slack 
bottom cable lengths cause only minor differences in the overall impact behavior and 
mechanical response of the fence. For example in the simulations with slack cables short 
duration force spikes in the fence impact force time histories were noted. 


6.3 Observations 


From the computational results of the analyses which are summarized in Table 2 and 3, the 
following general observations can be made: (i) for a fence with a bottom cable: elastic 
cable forces tend to reach a maximum when the impact takes place near, and above, the 
attenuator centroid, and (ii) for a fence without a bottom cable the elastic cable forces 
tend to increase with impact height. In both fences, with or without bottom cables it is 
observed that the rock post-impact velocity decreases with increasing impact height. 


The time histories presented in Figs. 14 and 15 give more detailed insight to cable 
behavior and rock velocity characteristics for a particular impact. Relevant observations 
are summarized as follows: 


(i) There exists a significant difference in the fence cable forces computed, assuming 
elastic and perfect elastic-plastic cable behavior, respectively. It is noted that the 
elastic-plastic material model predicts minor yielding which "cushions the impact", 
and reduces the unrealistically high impact peak forces in the cables obtained by the 
elastic cable model (see Figs. 14b, and 14f). 


(ii) The slack in the bottom cable produces force spikes introducing further impacts to the 
attenuators. Slack in the bottom cable causes a delayed sequence of impulsive cable 
forces which act between adjacent attenuators This phenomenon leads to force 
spikes in the cable forces (see Figs. 14c, and 14g). 


(iii) The fence with a bottom cable produces multiple impacts between the attenuator and 
the rock, thus, reducing the rock velocity gradually in stages. 


(iv) Minor plastification takes place in the top overhead cable mostly during the first 
impact between the rock and the attenuator (see Figs. 14i and 15i). 


7.0 CONCLUDING REMARKS 


A powerful numerical modeling tool has been developed for the non-linear dynamic 
impact modeling of a flexible rockfall prevention fence. An approximate three-dimensional 
impact analysis has been performed in an effective manner with a new two-dimensional 
DEM model. In this model the fence and rock kinematics are assumed to be planar, 
whereas the cable stiffness effects are modeled in a three-dimensional manner. Nonlinear 
cable effects have also been included. 


The DEM simulation procedure has been mathematically validated and applied to a series 
of full scale fence-rock impacts with varying design parameters. The numerical results for 
this series of impact simulations predict the maximum forces generated within the fence 
cables during a rock impact event. Corresponding predictions of post impact rock 
velocities have also been made. These initial calculations presented here illustrate the ease 
with which a designer can evaluate the effects of varying different fence design parameters 
on the fence capacity. It is hoped that further use of this numerical model, in conjunction 


with additional experimental data, will enable the design engineer to improve upon the 
current methodology employed in the design of flexible rockfall prevention fences. 
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Top Cable Span: 

Top Cable Sag: 

Modulus of Elasticity: 
Top Cable Area: 
Bottom Cable Area: 
Tire Height: 
Attenuator Height: 
Attenuator Width: 
Attenuator Unit Weight: 
Attachment Length: 
Rock Unit Weight: 
Rock Weight: 


Rock Impact Velocity: 
Coefficient of Restitution: 


L = 55 ft 

f= 5.5 ft 

E = 20,000 ksi 

Arp = 1-076 in? 

Ayot = 0-25 in? (doubly wrapped) 
t = 8 in 

12t = 8 ft 

w= 2.5 ft 

Yatr = 12.2 lbs/ft3 

d = 2.5 ft 


Yrock = 165 lbs/ft3 

600 lbs (case (i)) 

1200 lbs (case (ii)) 

Vy = 65 ft/s (horizontal) 
e= 0.4 


Table 1. Model Dimensions and Properties 
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FORCE IN 
FORCE IN TOP CABLE FORCE IN 


TOP JOINT WITH BOTTOM BOTTOM CABLE 
104 [Ibs] CABLE NO SLACK 105 [Ibs] 
105 [Ibs] 


WITH BOTTOM 


G. 83) (2: 00) (1. 5) 


Table 2. Maximum Forces for a 600 lb Rock Impacting a Fence with: (a) with Bottom 
Cable 0.5 ft. Slack, and (b) without a Bottom Cable 
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FORCE IN 
FORCE IN TOP CABLE FORCE IN 
TOP JOINT WITH BOTTOM BOTTOM CABLE 


104 [Ibs] CABLE NO SLACK 10> [Ibs] 
105 [ibs] 


WITH BOTTOM 
CABLE,SLACK = 0.50 FT. 


i 
(6.4) (2.00) 
(6.28) (2.00) 


(-17.8) 


Table 3. Maximum Forces for a 1200 Ib Rock Impacting a Fence with: (a) with Bottom 
Cable 0.5 ft. Slack, and (b) without a Bottom Cable 
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FORCE IN 
FORCE IN TOP CABLE FORCE IN Vx, - ROCK 
TOP JOINT WITH BOTTOM BOTTOM CABLE AFTER IMPACT 


104 [Ibs] CABLE NO SLACK 10> [Ibs] 
105 [tbs] 


WITH BOTTOM 
CABLE SLACK = 0.25 FT. 


-3.0 . . . . 
(5.77) (2.00) (1.00) (-2.3) 


(2. 00) a 00) (-7.2) 
(5. 89) (2. 00) (1. 00) (-3.7) 
(2. 00) (0. 89) (5.5) 


Table 4. Maximum Forces for a 600 lb Rock Impacting a Fence with: (a) with Bottom 
Cable 0.0 ft. Slack, and (b) (a) with Bottom Cable 0.25 ft. Slack. 
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FORCE IN 
FORCE IN TOP CABLE FORCE IN 


TOP JOINT WITH BOTTOM BOTTOM CABLE 
104 [Ibs] CABLE NO SLACK 1 
10> {Ibs] 


(4.84) (1.87) (1.00) (-8.4) 

5.82) (2.00 (1.00) (6.7) 

ik a 
(3.87) (1.68) (1.00) (-4.2) 
(6.68 (2.00) (1.00) (1.6) 
(6.33) (2.00) (1.00) (6.4) 
WITH BOTTOM 
CABLE SLACK = 0.25 FT. 


Table 5. Maximum Forces for a 1200 lb Rock Impacting a Fence with: (a) with Bottom 
Cable 0.0 ft. Slack, and (b) (a) with Bottom Cable 0.25 ft. Slack. 
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Fence structure - Three-Dimensional View 


Fig. 1 


26 


ob & ly 
ne i 


Top Cable 
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Bottom Cable 


(b) 


(b) Side View of a Single Attenuator and Rock 


Fence Structure: (a) Front View, 


Fig. 2 
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Attenuator Numbering System 


Fig. 3 


28 


Fig. 4 Attenuator Geometry 


Fig. 5 Cable Geometry - ith Cable Segment 
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Fig. 6 Contact Model: (a) "Normal" Contact Element, 
(b) “Shear™ Contact Element 
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Fig. 7 Particle Impact with Elastic Cables: 
(a) Velocities before Impact, 
(b) Velocities Immediately after Impact, 
(c) Post-Impact Velocities 
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DYNAMIC IMPACT GF TWO BODIES. e = @.5 
ROCK IMPACTING BODY WITH CABLE 


X - VELOCITY (M/S.) 


TIME IN (SECS } 


Fig. 8 Particle Impact with Elastic Cables: 
Velocity Time History 
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(a) (b) (c) 


Fig. 9 Pivoted Cylinder Impact Problem: 
(a) Velocities before Impact, 


(b) Velocities Immediately after Impact, 
(c) Post-Impact Velocities 
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SINGLE ATTENUATOR: FIXED PIVOT 
FENCE 48@ LA. ROCK 12@@ LBS. e = @.5 


TIME = -.BQ@E-B3 SEC. 


TIME = -3128 SEC. 


TIME - -6328 SEC. 


Fig. 10 Pivoted Cylinder Impact Problen: 
Geometry Plot 
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SINGLE ATTENUATOR: FIXED PIVOT 
FENCE 488 LB. ROCK 1208 LBS. e = @.5 
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Fig. 11 Pivoted Cylinder Impact Problem: 
Velocity Time History Plot 
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ROCK WEIGHT = 622 LBS 


TIME = +. 1B71E-@i SEC. 


-B4BBE-@1 SEC. 


od 
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Per 
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TIME = - 1888 SEC. 


Fig. 12 Impact of 600 LB Rock, with Bottom Cable (Slack = 0.5 FT): 
Geometry Plot 
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ROCK WEIGHT = 686 LBS 


TINE = +.1B71E-@1 SEC. 


-BABGE-B1 SEC. 


Fig. 13 Impact of 600 LB Rock, without Bottom Cable: 
Geometry Plot 
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Fig. 14 Impact of 600 LB Rock, with Bottom Cable (Slack = 0.5 FY): 
(a) to (d) Time History Plots for Elastic Cables, 
(e) to (i) Time History Plots for Elaso-Plastic Cables 
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Fig. 15 Impact of 1200 LB Rock, with Bottom Cable (Slack = 0.5 FT): 
(a) to (d) Time History Plots for Elastic Cables, 
(e) to (i) Time History Plots for Elasto-Plastic Cables 
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APPENDIX I - Modeling of Single and Double Layered Flexible Rockfall 
Prevention Fences 


The following technical paper describes the additional work required to modify the 
FENCE computer model for the analysis of flexible rockfall prevention fences 
constructed with two layers of attenuator columns. This publication also illustrates the 
capability of the FENCE model to simulate multiple rocks impacting a fence structure. 


DYNAMIC ANALYSIS OF FLEXIBLE STRUCTURES SUBJECT TO 
MULTIPLE ROCK-RUBBLE IMPACTS 


by 


G.G.W. Mustoe 
Division of Engineering 
Colorado School of Mines 
Golden, CO 80401, USA. 


Accepted for the Proceedings of the 2nd International Conference on Discrete Element Methods 
at 
M.LT., Boston Massachusetts, March. 1993 


Abstract 


A discrete element model for simulating the dynamic impact response of a highway rockfall prevention 
barrier system is described. This work is an extension of previous work for modeling single layer fence 
barriers subject to single rock impacts. The DEM technique employs two-dimensional modeling 
techniques to develop an efficient algorithm for a three-dimensional nonlinear dynamic structural system. 
The approach developed herein is suitable for modeling multi-layered flexible barrier systems impacted by 
single or multiple rocks. The resulting software can be used on microcomputers and is suitable for 
sensitivity studies and design calculations. 


INTRODUCTION 


Rockfall prevention highway barrier systems consisting of flexible fence structures are increasingly being 
used to minimize the risk of rockfall related automobile accidents or road closures. These barriers are 
placed near highways in mountainous regions wherever there is a significant potential rockfall hazard. 


Traditionally, most rockfall prevention highway barrier systems are rigid or semi rigid structures such as 
fences or walls. Rigid barrier systems are usually constructed from a combination of wood, metal grids 
and low cost geo-materials which are often available onsite. Because these barriers are designed to 
completely stop falling rocks they are usually fairly massive structures requiring a large amount of 
material and significant construction. Alternative rockfall prevention barrier systems are catchment pits or 
trenches which must be excavated or blasted so that they are deep enough to catch the rocks. 


More recently, highway engineers have developed various flexible barrier systems which overcome many 
of the problems associated with rigid structural barriers or pits. These flexible systems are much smaller 
and lightweight in construction, and are designed so that minor damage to the barrier can easily be 
repaired. It should be noted that minor damage should not significantly reduce the effectiveness of the 
barrier. The basic design concept of these flexible barrier systems is based upon the notion that the 
flexible barrier does not have to completely stop the falling rock, but just reduce the kinetic energy of the 
falling rock to an acceptable level after impacting the barrier. The kinetic energy of the rock after 
impacting the barrier must be such that the rock comes to rest harmlessly in a rolling zone. The concept of 
reducing the falling rocks kinetic energy is usually termed the principle of rockfall attenuation. 


In a flexible barrier system the rock velocity is slowed down by transferring some of its kinetic energy into 
the flexible barrier system in the form of impact energy losses, inelastic deformational energy and an 
increase of kinetic energy within the flexible barrier. Typical design considerations for a flexible barrier 
system are: (i) the accurate prediction of maximum impact forces within the flexible barrier system, and 
(ii) an estimate of the corresponding kinetic energy loss of the falling rocks. 


1.2 


In order to address these design issues full-scale testing of prototype flexible fences and barrier systems 
has been performed by the Colorado Department of Highways (CDOH). This experimental work has been 
successful in determining the overall performance of typical flexible barrier systems, estimating the basic 
design parameters, and modeling barrier failure mechanisms. The major limitations of the data obtained 
from this test program are: (i) it is difficult to scale the experimental data for much larger or smaller 
structures, (ii) limited quantitative information of falling rock kinetic energy losses, and (iii) limited 
applicability to modifications in the barriers structural design, and (iv) most of the data is qualitative . 


Recently, a numerical model based upon a combination of the discrete element methodology and a 
conventional dynamic structural model has been developed to simulate the dynamic response of a 
particular flexible rockfall prevention barrier system [1]. This barrier type consists of hanging columnar 
attenuator masses assembled from used truck tires and rims, threaded onto a steel rod, which are 
suspended from a horizontal overhead wire cable and connected by a bottom cable. Fig. 1 (a) illustrates a 
typical two fence layer flexible barrier system. The bottom cable connection is installed to increase the 
kinetic energy losses of the impacting rocks. This increase occurs because the combined mass of all 
attenuators are more tightly linked together. It should be noted that any slack in the bottom cable results 
in a staggered motion of the attenuators. When a rock hits the barrier the hanging columns begm to move 
in 2 pendulum-like motion. The sizes of the most likely impacting rocks range between 500 Ibs to 1500 
Ibs. Corresponding impact rock velocities vary between 50 and 70 ft/s. It should be noted that larger 
impacting rocks that are over 2000 Ibs in weight are likely to cause significant fence damage, and exceed 
the barrier system capacity. 


The numerical model described here is a technique that can be used in conjunction with experimental 
methods to estimate the barrier capacity. For example, the dynamic impact calculations can study the 
effects of different rock sizes and rock velocities on particular barrier configurations. Barrier response 
parameters of particular importance include; (i) the forces in the hanger connection and the overhead and 
bottom cables, and (ii) the decrease in rock velocity due to impact. 


The numerical procedure developed is based upon the DEM. This procedure solves the nonlinear coupled 
dynamic equations of motion for the idealized barrier structure and the impacting rock. The DEM model 
assumes that the barrier structure can be idealized as a system of connected rigid bodies, and the rock as 
an impacting rigid body. The dynamic barrier-rock impact conditions are modeled with a penalty 
formulation and an automatic contact detection algorithm. This numerical procedure utilizes an explicit 
time stepping scheme to solve the discretized equations of motion, see [1] for further details. 


The work presented in this paper is an extension of the analysis technique and associated computer 
software developed by Mustoe and Huttlemaier [1]. The new developments reported here include: (i) the 
modeling of rockfall events consisting of multiple rock impacts, and (ii) dynamic analysis of flexible 
barrier systems consisting of multiple fence structures. 


BARRIER DEM MODEL 
The rockfal! barrier discrete element algorithm includes the following features: 


(® an automatic contact detection scheme which automatically recognizes new contacts as the simulation 
progresses, and updates the contact topology appropriately, 

(ii) a contact force model that determines the forces acting between interacting fence attenuator columns, 
impacting rocks and connecting fence cables, and 

(iii) a time integration scheme to solve the governing equations of motion of the fence attenuator columns 
and impacting rocks, that allows for finite displacement and rotation. 


1.3 


The modeling aspects and assumptions of the discrete element model, namely the rigid body dynamics of 
the barrier system and rock, the cable behavior, the rock/attenuator impact interaction model and the 
numerical integration of the dynamic equations of motion are discussed below: 


The impacting rocks and fence attenuator columns are modeled as three-dimensional rigid spherical 
and cylindrical bodies respectively, that are confined to move in a planar manner within the x,y plane 
(see Fig 1 (b)). 

A superquadric representation [2] of the two-dimensional boundaries of the impacting rocks and 
fence attenuator columns is employed. This kind of geometrical modeling has been used previously by 
Williams and Pentland [3] in DEM analysis. 

Multiple fence layers within the flexible barrier system interact with each other because of the contact 
between adjacent attenuator columns in different fence layers. This phenomena occurs because of the 
planar kinematics of the attenuator columns. 

Multiple interacting rocks are assumed to hit the central localized portion of the barrier that is limited 
to one or two attenuators in the center of the barrier. This assumption is justified since the diameter of 
most impacting rocks is less the diameter of the cylindrical attenuator columns. _ 

The connecting cable segments between neighboring attenuators in the same fence layer are modeled 
as elastic or elasto-plastic massless elements that only transmit tensile forces. 

Contact forces are generated from contact stiffness and damping elements which act between the 
central attenuator(s) and the impacting rocks or between adjacent attenuator columns in different 
fence layers. These forces are determined from the current penetration or overlap distance and the 
relative velocities of the bodies in contact at the common point of contact. Note, that contact forces 
between adjacent attenuator columns in the same fence layer are neglected, since they are small 
compared with the impact force between the central attenuator and the impacting rocks. 


BARRIER MODEL EQUATIONS 


Equations of Motion 


The equations of motion for the ith fence attenuator within a fence layer modeled as a rigid body, 
assuming planar dynamic motion confined in the x-y plane are given by: 


m; Xx} = ZF, (1) 
my; = ZF; (2) 
1,6, ==M, (3) 


for i = 1,2,...N, where N, is the total number of fence attenuators, and N is defined as N= WN, + 1/2. 
This definition of N is because of the impacting rock(s) are assumed to hit the central fence attenautor 
column. See Fig. 1 (c) for details of the attenuator DEM idealization. The corresponding resultant loading 
terms in equations (1),(2) and (3), are: 


LF, = Ty + By + P, (4) 


1.4 


ZIM, = B, beos0; + By, bsin®; + T,,acos@; + T, asind; 
+P, Ay; + Pyi Ax; 


(6) 


where T,; and Ty are the global x- and y- components of tension forces induced at the top overhead cable 
connection, B,; and By are the global x- and y- components of tension forces induced at the bottom 


cable connection, P,; and P,; are the global x- and y- components of the contact forces generated by the 
rock attenuator impact, and g is the gravitational constant. For the ith attenuator the following parameters 
are defined: m, is the mass, I, is the mass polar moment of inertia with respect to the centroid, x; and y; 
are the global centroidal coordinates, and 0; is an angular coordinate defining the axial direction of the 
attenuator. The distances of the top and bottom cable conzections with respect to the attenuator centroid 
are defined by a and b, respectively. The global x- and y- coordinates of the point of contact of the rock 
with the attenuator measured with respect to the attenuator centroid are defined by Ax; and Ay;. Note, 
that if N is even, the mass and moment of inertia terms in equations (1) through (3), for 1, are doubled. 
This modification is required because of the assumed symmetrical impact geometry which implies that the 
two central attenuators are hit simultaneously. 


The corresponding equations of motion for impacting rocks are similar to those developed above for a 
fence attenuator, and are not repeated for sake of brevity. Note, that the cable force terms in the equations 
of motion for the rocks are zero. 


Cable Attachment Forces 


The top and bottom cable attachment segment forces are computed using simple equivalent axial cable 
stiffness coefficients and determining the current change in cable segment lengths from their original 
initial segment lengths. It is important to note that the initial tension in the cables must be accounted for 
in these force calculations. For further details, see Mustoe and Huttlemaier [1]. When modeling the slack 
in the bottom cables it should be noted that a significant amount of slack cable length must be taken up 
before any tensile forces are transmitted to the attenuators. 


Interaction Force Laws 


A simple dynamic contact model is used, since the energy transfer mechanism between the rock and the 
fence is of primary interest. The normal contact element model is comprised of a linear spring stiffness 
and viscous damper, connected in parallel, that acts in the contact normal direction at the instantaneous 
point of contact. Behavior in the contact shear direction is defined by a shear contact element model 
consisting of a linear spring stiffness and a slider that simulates a simple Coulomb dynamic friction 
model. Note, that the contact normal and shear directions are defined by the boundary geometry of the 
body with the smoothest tangent plane at the point of contact, The normal and shear contact element 
model are active only when the discretized boundaries of the rock and the attenuator penetrate each other. 


The normal instantaneous contact force, F_, is defined by: 


F, = 


k, 5, +08, 6, 20 
0 6. <0 


r.5 


where, k,, is the normal contact stiffness, c is the normal contact viscous damping coefficient, and 5, is 
the normal penetration distance at the point of contact. 


The shear instantaneous contact force, F,, is given by: 


k,5, —_k,|8,| <w|F,] and 8, >0 
F=(uF,  &k,|6,|2u/F| and 6, >0 (8) 
0 §, <0 


where, k, is the shear contact stiffness, 5, is the shear slip distance at the point of contact and ji is the 
dynamic coefficient of friction between the rock and the attenuator. 


Note: A suitable value of the normal contact stiffness k, is specified by the user by limiting the 
maximum penetration, or overlap, between the surfaces of the impacting two bodies to a predefined value. 
Typically, this maximum penetration is defined as 2 small fraction of the average body dimension of the 
contacting bodies. An estimate of k, cam simply be determined with an approximate energy calculation 
for the two impacting bodies during an elastic collision. This approximate calculation assumes that the 
bodies are approximated as two particles of mass m, and m,, respectively. Using the above approach an 
estimate for k,, is given by: 


k, = [ m, my /(m, +m) Wie /Avaux (9) 


where v,,, and A, are the estimated maximum normal relative speed and overlap distance at the 
contact point respectively during the impact. 


Restitution Model 


In the discrete element numerical model, the energy loss due to impact is simulated with a viscous damper 
which acts at the point of impact, P, in the contact normal direction. The normal viscous force in this 
damper, which is the second term in the normal contact force equation (7), is given by: 


FE, =c8 (10) 


where é, is the normal component of the relative velocity of the fence with respect to the rock at the 
point of contact, and is defined by: 


5, =(F_—- ¥,,) fi (il) 


Note, that fi is the unit outward contact normal direction at the point of impact, P, and ¥p, and ¥p, are 
the velocities of the fence and the rock during the impact, respectively, and c is the normal contact viscous 
damping coefficient as defined in equation (7). 
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In a dynamic impact between two bodies, the energy loss can be approximately modeled with a simple 
restitution model that requires e, the coefficient of restitution, to be known apriori via experimental data. 
Due to lack of available quantitative data concerning the impact energy losses during a barrier-rock 
impact situation this simple model is appropriate. However, in order to apply the normal contact force 
law, equation (7), a relationship between the viscous damping coefficient c, and the coefficient of 
restitution ¢, is required. For this purpose, a simple approximate relationship between the viscous 
damping coefficient, c, and the coefficient of restitution, e, can be derived by assuming that the rock and 
the attenuator are approximated as two particles of mass m, and m,, respectively. This results in the 
relationship: 


c= 2 In(i/e)[k, m, mr /(m, +m,)] a [x? + (in(1/e))?]” (12) 


where k, is the normal contact stiffness at P, the point of impact. 


It should be noted that the use of equation (12) is only approximate when applied to a non-central impact 
between a spherical rock and a cylindrical attenuator. Some numerical experimentation, therefore, is 
usually required to obtain the precise value of the viscous damping coefficient which corresponds to a 
specified value of the coefficient of restitution. 


Contact Detection Algorithm 


The multiple contacts which occur between the rocks and the fence attenuator columns are detected 
automatically with an algorithm which performs the following functions: 

(i) At the beginning of the analysis the problem space is subdivided into a grid of rectangular cells or 
subregions. The dimensions of these cells is defined to be a function of the average dimension of the 
bodies within the problem and selected to minimize computational effort. 

Gi) Simple geometrical checks are carried out to determine the cell locations for each body in the problem 
throughout the simulation as the bodies move. 

(iii) At every time step local body-to-body contact checking is performed on a cell by cell manner for cells 
within the problem space. This algorithm is similar to that developed by Preece and Taylor [4], and the 
method employed by Zhang [5] for modeling particle flow problems. 


Explicit Time Stepping Scheme 


The time stepping scheme employed in the DEM procedure is the standard explicit central difference 
method [1}, applied to the equations of motions defined by equations (1) through (3). This leads to three 
sets of generalized velocity and displacement update formulas for the centroidal coordinates of a body. 
Large motion of the bodies ( the fence attenuators and rocks) is simulated with successive incremental 
updating of the bodies current coordinates. 


DEM IMPACT ANALYSES 

Model Dimensions and Properties 

A full size barrier model is used to study the dynamic response of the barrier system due to impacting 
rocks. A similar full size prototype barrier system has been built by the Colorado Department of Highways 
(CDOH). Limited experimental data, which is mainly qualitative, has been obtained from this barrier 
structure by the CDOH. The dimensions and material properties for these numerical analyses are as 
follows: . 


I.7 


Span of top cable 55 ft 

Sag of top cable 5.5 ft 

Youngs modulus of cable 20,000 ksi 
Cross section area of top cable 1.076 in2 
Cross section area of bottom cable 0.25 in? (double wrapped) 
Fence attenuator height 8 ft 

Fence attenuator width 2.5 fi 

Rod attachment length 2.5 fi 
Fence attenuator unit weight 12.2 Ibs/ft3 
Rock diameter 2.5 ft 
Rock weight : 1200 Ibs 
Rock impact velocity 53 fi/s2 
Coefficient of restitution _ 0.4 
Maximum estimated contact overlap 0.5 ft 


In the present computer model the cable behavior is modeled with a simple linear elastic model. It should 
be noted that more realistic material models for the cable can be readily inserted into the present 
algorithm. For example, a constitutive model with time dependent stiffness properties which accounts for 
cable life [6], could be easily implemented. 


Analysis Results 


A series of analyses were performed to study the influence of the following design variables: (i) barriers 
systems consisting of one and two fence layers, (ii) single and multiple rocks impacting the barrier 
systems, and (iii) barriers with and without a bottom cable. All analyses assumed that the bottom cable 
segments had a 0.5 ft initial slack length. The maximum estimated contact overlap of 0.5 ft was chosen to 
approximate the local contact deformation of the fence attenators which are constructed from rubber tires 
and rims. 


Each computer run was performed with a time step length of approximately At =2 x 10~ seconds. The 
analyses performed here required between 5 and 15 minutes of CPU, on an Inte! 80386, 33 MHz CPU 
microcomputer. 


Geometrical snapshots of the barrier geometry and impacting rocks for four different simulations are 
shown in Figs. 2 and 3. Figs. 2 (a) through (c) illustrate a single rock impacting a single layer fence 
barrier system with top and bottom cable attachments. Figs. 2 (d) through (f) illustrate a single rock 
impacting a two layered fence barrier system with top and bottom cable attachments. Figs. 3 (a) through 
(c) illustrate nine rocks (with a total combined weight of 1200 Ibs) impacting a two layered fence barrier 
system with top and bottom cable attachments, Figs. 3 (d) through (f) illustrate a single rock impacting a 
two layered fence barrier system with only a top cable attachment. 


Comparison of Fig. 2 (a) through (c) and Fig. 2 (d) through (f) shows the following: (i) similar dynamic 
rock motion during these two impacts, and (ii) the dynamic response of the two layered fence bartier is 
significantly less than that of the single fence barrier system. The simulations for the two layered fence 
barrier systems with one and two cable attachments respectively that were subject to a single rock impact 
illustrate the effect of the bottom cable attachment. The bottom cable stiffens the barrier system and 
decreases the rock velocity more than a barrier system without a bottom cable (see Figs. 2 (d) through (e), 
and Figs. 3 (d) through (e)). The effect of multiple or single rocks impacting a barrier system is illustrated 
in Figs. 2 (d) through (f), and Figs. 3 (a) through (c). The overall response of the barrier system is similar 
in these two cases. 


A comparison of the rock velocities for the above simulations are illustrated in Figure 4. All of the barrier 
systems with top and bottom cable attachments essentially absorb all the kinetic energy of the impacting 
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tock (or rocks). The barrier system with only a top cable only reduces the rock velocity from 53 ft/s. to 32 
ft/s. Figs. 5 (a) and (b) show the time histories of the dynamic forces generated in the top joint of the 
central vertical attenuator rod for the four simulations. Fig. 5 (a) shows the effect of the bottom cable in a 
barrier system on the dynamic loads generated in the top joint. The high peak load at 0.05 seconds in the 
barrier with two cables is because the 0.5 ft.of initial slack in the bottom cable has been taken up and 
become taut. Fig. 5 (b) shows that the highest peak force is experienced by the single fence barrier system. 
This figure also shows that the maximum force generated in the two layered fence barrier impacted by 
nine rocks is similar in magnitude to the single rock impact simulation. Note the effect of the multiple 
rock impacts is seen in the smaller peak loads after 0.08 seconds. 


CONCLUDING REMARKS 


A powerful numerical modeling tool has been developed for the non-linear dynamic impact modeling of a 
flexible rockfall prevention barrier system. An approximate three-dimensional impact analysis has been 
performed in an effective manner with 2 new two-dimensional DEM model. In this model the barrier 
system and rock kinematics are assumed to be planar, whereas the cable stiffness effects are modeled in a 
three-dimensional manner. Nonlinear cable effects can be modeled and have been discussed elsewhere [1]. 


A DEM simulation procedure for modeling flexible barrier systems with multi-layered fence structures 
subject to single and multiple rock impacts has been developed and demonstrated. This numerical model 
has been applied to a series of full scale barrier-rock impacts. The numezical results for this series of 
impact simulations predict the maximum forces generated within the fence cables and the column 
connections during a rock impact event. Corresponding predictions of post impact rock velocities have 
also been made. These initial calculations presented here illustrate the ease with which a designer can 
evaluate the effects of varying different barrier system design parameters on the barrier capacity. It is 
hoped that further use of this numerical model by the CDOH, in conjunction with additional experimental 
data, will enable the design engineer to improve upon the current methodology employed in the design of 
flexible rockfall prevention barrier systems. 
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Fig. 1 A Typical Flexible Rockfall Barrier System. (a) Three-dimensional View of a Two Layered Fence 
Barmier, (b) Side View of a Two Layered Fence Barrier, (c) DEM Idealization of a Fence Attenuator 
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SINGLE ROCK IMPACT.TWO CABLES 
ONE FENCE. 1288 LBS ROCK 
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Fig. 2 Geometrical Snapshots of a Barrier-Rock Impact Simulation, Figs. 2 (a),(b) and (¢) are for a Single 
Layered Fence Barrier, Figs. 2 (d),(e) and (f) are for a Two Layered Fence Barrier 
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Fig. 3 Geometrical Snapshots ofa Barrier-Rock Impact Simulation, Figs. 3 (a),(b) and (c) are for a Two 
Layered Fence Barrier subject to a Multiple Rock Impact, Figs. 3 (d),(e) and (f) are for a Two Layered 
Fence Barrier with a Top Cable Attachment only. 
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Fig. 4 Rock Velocity Time Histories for Different Barrier-Rock Impact Simulations 
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Fig. 5 Dynamic Forces generated in the Top Joint of the Central Vertical Attenuator Rod. (a) Simulations 
of One and Two Layered Fence Barrier Systems with Bottom Cable Attachments, (b) Simulations of Two 
Layered Fence Barrier Systems subject to a Multiple Rock Impact, and with a Top Cable Attachment only. 
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Computer Instructions and Users Manual for FENCE 
Abstract 


This manual is a detailed guide for the FENCE discrete element computer 
software package. It also describes the necessary steps that an engineer 
must undertake to perform a dynamic impact rock-fence analysis using the 
FENCE discrete element computer software package. 
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1.5.3 Sample Geometry Plotting Session 
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I.1 Introduction 


FENCE is a special purpose computer program based on the Discrete Element Method 
(DEM) written for the Colorado Department of Highways (CDOH). The program 
performs rock impact calculations for a specific type of flexible rockfall fence type, The 
fence model used for this analysis of a series of attenuators suspended from a steel cable 
and is based on a prototype developed by CDOH. The theoretical foundation of the 
discrete element formulation for this specific application is discussed in previous chapters 
of this report. The fence model description is also included earlier in this report, in 
addition, fence input parameters are discussed below. This Appendix is intended as a 
guideline for the program user and includes: (i) an overview of the computer software 
including the different executable modules, and the associated input, output and graphics 
data files, (ii) the general strategy required to perform an impact analysis using FENCE, 
(iii) a detailed description of the parameters within a FENCE input file, and (iv) an 
example that guides the user through the analysis procedure. 


IL.2 Program Overview and File Definitions 


The DEM dynamic fence-rock analysis software FENCE consists of three major modules: 
(i) the computational DEM executable module (FENCE) that performs the dynamic 
impact calculation, (ii) the interactive graphics geometrical snapshot and animation 
executable module (GPLOT), and (iii) the interactive graphics executable module for time 
history results (HPLOT). 


For the DEM computational module, FENCE, the input data contains the following 
parameters: 


(i) fence and rock geometry and mass properties, 

(ii) material model and properties for the fence cables, 

(iii) impact and restitution parameters between the fence and the rock, 

(iv) the initial rock velocity prior to impact, and 

(v) control parameters which describe the type and time interval sampling interval for 
saving computed results. 


This input data is contained in a data file with a user-specified name appended with the 
extension .DAT. 


The computed results are written to three different output data files with a user-specified 
name appended with the extensions OUT, .HST, and .GEO, respectively. These output 
data files are called the printout file, the time history file, and the geometry file, 
respectively. 


The printout file is used directly to verify the input data and for checking purposes by the 
user. The time history and the geometry results data files are input data for the graphical 
post-processing modules HPLOT and GPLOT, respectively. 


The HPLOT graphics module interactively displays time history plots and produces 
postscript graphics output. Typical time history data plots include; (i) centroidal 
Kinematics quantities such as position and velocity components for the rock and the 
attenuator columns, and (ii) forces in the overhead and bottom cable segments, and the 
attenuator hanger connections. 


The GPLOT graphics module is used to display the following geometrical data for the 
fence-rock system; (i) single snapshots of the fence-rock geometry at prescribed times 
during the impact, or (ii) animation of the fence-rock motion throughout the impact. 
Postscript graphics output of the single geometry snapshots can also be produced by 
GPLOT. 


In the example illustrated in Figure II.1 the user-specified name for the input data file is 
fname, and the user-specified name for the output data files is fname2. 


fname2.OUT 
(printed output) 


fname2.HST 
(time history output) 


fname2.GEO 
(geometrical output) 
(interactive graphics and hardcopy) (interactive graphics and hardcopy) 


(a) Dynamic Impact Calculation (b) time history graphics postprocessing (c) geometry snapshot postprocessing 


fname1.DAT 
(input data) 


Figure II.1 Execution Sequence of FENCE Computer Software: 


The data file definitions for the above figure are as follows: 


(i) fname] .DAT: Input File, where fname] is a user defined name, 

(it) fname2.OUT: Output File in ASCII format, where fname? is a user defined name, 
(iti) fname2.HST: Time History Plot File in binary format, 

(iv) fhame2.GEO: Geometry Plot File in binary format, 

(v) hplotfile: 'Postscript' File for printer plot of time history results, 

and (vi) gplotfile: ‘Postscript' File for printer plot of geometry snapshot results. 


IL3 Modeling Strategy with FENCE Software 


I1.4 


When performing a rock-fence impact analysis the following strategy is recommended: 


(i) Gather the relevant design data about the fence structure and the impact conditions. 
This data includes: (a) the geometrical configuration of the fence and attenuator masses, 
and (b) the shape and mass of the impacting rock and its initial impact position relative to 
the fence and the impact velocity. 


(ti) Make a simple sketch of the side view ( x,y plane) of the impact geometry. 
This will be used to define the input geometry data. 


(iii) Decide what results are required from the analysis. This includes the following: (a) 
how much printed output is needed ( This can be voluminous, so be careful.), (b) which 
elements (fence, rock and attenuators) must be flagged for results output purposes, (c) 
what time history information is needed, and (d) what geometrical snapshot information is 
needed. 


(iii) Generate an input data file with a user defined name with an extension .DAT. 


(iv) Perform an initial test run to verify the input data. This can accomplished by running 
the analysis for one time step only. Check the input data by reviewing the .OUT file for the 
following likely errors: (a) incorrect geometry of the fence, or initial position of the 
impacting rock, (b) check fence and rock masses and initial rock velocity, (c) check the 
total span and initial sag of the top overhead cable, (d) check the horizontal spacing 
between attenuator columns, (e) check maximum specified overlaps between rocks and 
attenuators elements and restitution coefficient. Verify the problem geometry by viewing 
the initial configuration with the GPLOT graphics module. 


(v) Re-run problem for approximately 200 time steps to check the initial impact 
kinematics, such as the initial speed and direction of the impacting rock. Check the initial 
impact parameters graphically with GPLOT. If there is a problem with the impact 
parameters change the initial impact parameters in the .DAT input file. If the data verifies 
determine the number of time steps required for a complete simulation, change the DAT 
input file, and make the final nun for the complete impact simulation. 


NOTE: There are number of example data files provided with the DEM computer 


software. Therefore, many of the above steps can be performed quickly by copying an 
existing data file and modifying it accordingly for a particular rock-fence impact analysis. 


11.3.1 Details of the Fence-Rock Geometrical Configuration Definitions required by 
FENCE Input Data 


In a typical fence-rock impact analysis only the geometry of three bodies is required. 
These bodies (or elements) are usually: (i) the ground or base element, (ii) the rock 
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element, and (iii) the central attenuator element. Figure II.2. shows the general situation 
for a single rock impacting a fence. For the problem defined in Figure 1.2 it should be 
noted that the geometry for elements 1 through 3 has to be defined. The geometry of 
elements 4 through 7 is generated automatically. Note that elements 4 through 6 are slave 
attenuator elements, and element 7 is a dummy fixed attenuator element . 
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(a) X-Y Plane View (o) Front View 


Figure [1.2 Typical Geometrical Configuration for a Rock - Single Layer Fence Impact 
Analysis 


The following points should be noted in defining the geometry for a problem: 


(i) the elements representing the ground, the rock and the central attenuator should be 
defined in this order. For example, the ground is element 1, the rock is element 2 and the 
central attenuator is element number 3. Note that if multiple rocks are defined they must 
be defined sequentially. In a three rock impact analysis, the ground is element 1, the rocks 
are elements 2 through 4, and the central attenuator is element number 5. 


(ii) the total number of attenuators, not accounting for symmetry must be defined in the 
input data. For the problem defined in Figure 2, the total number of attenuators is 7. This 
number is computed from the central attenuator and three slave attenuators on both sides. 
Note an additional attenuator element is used for fixity purposes in the FENCE model. In 
this example it is element 7. This element, however, is not defined by the user. 


(iit) the total span, S, the sag height h, and the spacing between the centers of adjacent 
attenuators d, must be defined 


(iv) the ith cable segment is defined as the segment on the immediate right of attenuator 
element 1. 


(v) an example of the element numbering for a two layered fence barrier is shown in 
Figure 11.3. 
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(a) Front View (2nd Row) (b) X-Y Plane View (c) Front View (ist Row) 


Figure 11.3 Typical Geometrical Configuration for a Rock - Double Layer Fence Impact 
Analysis 


113.2 Boundary Geometry Definition of the Ground, Attenuators and Rocks with 
Superquadrics 


In the two-dimensional discrete element FENCE model the rigid bodies that represent the 
attenuators, the ground surface and the rock(s) are idealized with a superquadric 
geometrical description. For further details of superquadric functions, see Barr (1981). 


For the ith body the following parameters are defined: m, is the mass, I, is the mass polar 

moment of inertia with respect to the centroid, X, and Y, are the global centroidal 
coordinates, and @. is an angular coordinate defining the direction of the local x-direction 
defined by the principal axis of the body with respect to the global X-direction (see Figure 
11.4). The parameters that define the shape of the body are aj, b;, and n; and are defined as 
the semi-lengths of superquadric in the local centroidal x- and y- directions, and the real 
exponent which specifies the roundness of the superquadric respectively. 


They 


a ‘, 
Figure 1.4 Definition of Local and Global Superquadric Geometry 


IL4 Description of the FENCE Input Data File 
This section describes the input parameters and the required format for the input file with 


an extension .DAT. The READ statements described below are contained in the input 
source file INPUT.FOR. 


READ(*,'(A)') flinp 

flinp = screen prompt for input data file (e.g. fnamel) 
READ(*,'(A)’) flout 

flout = screen prompt for output files (e.g. fname2) 
READ (iin, '(a)') ptitle 

ptitle = title line for input data (1st line in fnamel.DAT) 


READ(in,*) fact_dt.niter.nprint.nhisto.f penet 


fact_dt = factor of critical time step length 

niter = number of time steps (iterations) 

nprint = number of time step intervals between writing results onto .OUT and 
.GEO files 

nhisto = number of time step intervals between writing results to HST file 

f_penet = minimum penetration between rock and attenuator during impact defined 


as fraction of the bounding dimensions of bodies. 


NOTE: This can be overidden by using an explicit definition of maximum penetration 
between attenuator and rock elements and attenuator and attenuator elements. See 
READ(iin,*) rf_overlap, .... 


READ(in,*) nel,igrav.grav 
nel ~=number of elements (e.g., for slope, attenuators and rock (or rocks)) 
NOTE: usually nel = 3 
igray = 1 ... gravity acts in positive y-direction 


igray =-1 ... gravity acts in negative y-direction 
grav = gravitational constant 


READGin,*) elstore 


elstore = number of elements tagged for recording time history results 


READ(iin,*) (elnos(i),i=1,elstore) 


elnos = array of element I.D's tagged for recording time history results 


READ(in,*) totvar 
totvar = no. of variables to be stored (max. =6) 
READ (iin,*) (variab(i),i=1 bvmax) 
variab = list of time history variable codes to be recorded for elements in elnos() 
NOTE these integer codes are defined in global._par source file and printed in the COUT 
e. 
*** beginning of element loop, i = 1,nel 


READ (in * sha oints 


type = (rock=0, attenuator=1 ) 

shape = (sphere or 
( ellipsoid (rock) =0) 
( rectang. block (rock) =1) 
(cylinder (attenuator) = 2) 


points = number of nodal points on superquadric 
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READ(iin,*) xlen,ylen,power,rho,barlen 


xlen = length in x-dir. (a; ) 

ylen = length in y-dir. (b;) 

power = power of superquad (nj) 

rho = mass density 

barlen = attachment length (attenuator only) 


Note: See Figure 11.4 for superquadric geometry definitions 


READ(in,*) xcent,ycent,theta omega.f: 
initial global geometrical and dynamic centroidal conditions for an element: 


xcent = x-coordinate 
ycent = y-coordinate 
theta = rotation 

vx = x-velocity 

vy = y-velocity 

omega = angular velocity 
fx = force in x-direction 
fy = force in x-direction 
mz = moment 


READ (in,*) be(1),be(2),.be(3),x1_fix.yl_fix 


be(i) = centroidal fixities 
be(i) = 1 .. fixed 

be(i) = 0 .. free 

be(1) = x-direction 

be(2) = y-direction 

be(3) = rotational 


*** end of element loop 


READ(in,*)span.sag.dcol,youn to bot,nocab 
ncol.ecoeff,area_top.area_bot,frict.slack,brek top, 


brek bot 


Definition of fence layer parameters: 
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span = Horizontal cable span 


sag = Vertical cable sag 

dcol = Horizontal distance between columns 

young = Youngs modulus of cable 

y_top = y-coordinate of top cable connection 

y_bot = y-coordinate of bottom cable connection 
nocab = No. of top/bottom connections on attenuator 
ncol = No. of tire attenuator columns in a fence layer 
ecoeff = coefficient of restitution 

area_top = x-section area of top cable 

area_bot = x-section area of bottom cable 

frict = coefficient of friction 

slack = slack in between bottom connections 
brek_top = top cable breaking strength 

brek_bot = bottom cable breaking strength 


READ(iin,*) no_fence, sp_fence 
Multiple fence layer parameters: 


no_fence = number of fence layers 
sp_fence = center to center spacing between the two fence layer 


READ(iin,*) rf_overlap, ff overlap 


Definition of maximum penetration between attenuator and rock elements and attenuator 
and attenuator elements in different fence layers: 


rf_overlap = approximate maximum penetration between attenuator and rock 
ff_overlap |= approximate maximum penetration between attenuator and attenuator in 
different layers 


*** End of Input Data file *** 


IL5 Example Aualysis Datafiles and Results 


IL5.1 Sample Input Datafile 
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The following data set is for a double layer fence structure being impacted by a single 
rock. This data set is called sample2.DAT on the diskette supplied with FENCE. 


Data Optional Comments 

DOUBLE LAYER FENCE ROCK WEIGHT = 1200 LBS 

0.0 3000 60 30 0.05 fact_dt,niter,nprint,nhisto,fact_penet 
3 -1 32.18 nel, igrav, grav 

5 elstore 

2,3,4,5,6 (elnos(i),i=1,elstore) 

8 totvar 

51 52 48 67 8 35 360 (variab(i),i=1,hvmax) 

0.0 1.0 2.0 type,shape,points - SLOPE ELEMENT 
10.0 0.5 0.0 0.10.5 xlen,ylen,power,rho, barlen 

0.0 -5.0 0.0 0.0 0.0 0.0 0. 0. 0. xc,yc,th, vx, vy,oma, fx fy,mz 

1.0 1.0 1.0 0.0. be(1),be(2),be(3), xfix, yfix 

0.0 0.0 4.0 type,shape,points - ROCK 

1.201 1.201 0.0 5.13 0.0 xlen, ylen, power,rho, barlen 

2.451 0.0 0.0 -65. 0.0 34. 0. 0. 0. x¢,yc,th, vx, vy,oma, fx, fy,mz 

0.0 0.00.0 0. 0. be(1),bc(2), bc(3), xfix, yfix 

1.0 2.0 2.0 type,shape,points - FENCE attenuator 
1.25 4.0 0.00.76 2.67 xlen, ylen,power,rho,barlen 

0.0 0.0 0.0 0.0 0.0 0.0 0. 0. 0. xc, yc, th, vx, vy,oma, fx, fy,mz 

0.0 0.0 0.0 0.0. be(1),be(2),be(3),xfix, yfix 

55.0 5.5 2.5 2.0e7 6.67 -4.0 2 13 0.4 1.076 0.25 0.0 0.5 le8 1e8 

2, 2.6 - no. fences, spacing between fences 
0.6, 0.5 - max. rock and fence overlaps 


span,sag,dcol,e,y_t,y_b,nocab,ncol,ecoeff,a t,a_b,frict,slack,brek_top,brek_bot 


***e% End of DATA ***** 


11.5.2 Sample Output Datafile 


The following is selected printed output from sample2.OUT that was produced from 
executing FENCE with the data set sample2.DAT. This output data set is on the diskette 
supplied with FENCE. 


DOUBLE LAYER FENCE ROCK WEIGHT = 1200 LBS 
factor of critical time step = .000E+00 

total no. of time steps = 3000 

time steps per printout = 60 

history time steps = 30 
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no. of discrete elements = 3 
igrav ( =1 for positive y-dir.) = -1 
gravity = .3218E+02 


no. of etements for time history records= 5 
element history array: 

23 4 5 6 
no. of variables to be stored (max. 8) = 8 


list of variables to be stored: 
(see pointers defined in global.par) 


51 52 48 6 7 8 35 36 


Properties for element no.: 1 
Type (rock=0,fence=l1) = O 
Shape = ] 

( sphere or ) 

(ellipsoid (rock) =0) 

( rectang. block (rock) = 1) 

(cylinder (fence) =2) 
Number of points on superquad = 2 
Length in x-dir. (large radius)= .100E+02 
Length in y-dir. (short radius)= .500E+00 
Power of superquad = .000E+00 
Mass density = .100E+00 
Attachment length (fence only) = .500E+00 


Properties for element no.: 2 
Type (rock=0, fence=1) = 0 
Shape = 0 

( sphere or ) 

(ellipsoid (rock) =0) 

( rectang. block (rock) = 1) 

(cylinder (fence) =2) 
Number of points on superquad = 4 
Length in x-dir. darge radius)= .120E+01 
Length in y-dir. (short radius)= .120E+01 
Power of superquad = .000E+00 
Mass density = 513E+01 
Attachment length (fence only) = .000E+00 


Properties for element no.: 3 
Type (rock=0,fence=1) = 1 
Shape = 2 
( sphere or ) 
(ellipsoid (rock) =0) 
( rectang. block (rock) = 1) 
( cylinder (fence) =2) 
Number of points on superquad = 2 
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Length in x-dir. (large radius)= .125E+01 
Length in y-dir. (short radius)= .400E+01 
Power of superquad = .000E+00 
Mass density = .760E+00 
Attachment length (fence only) = .267E+01 


Element initial conditions 

el. x-loc y-loc theta vx vy ang-v fx 
mz 

1 .000E+00 -.500E+01 .000E+00 .000E+00 .000E+00 .000E+00 .000E+00 .000E+00 .000E+00 
2 .245E+01 .000E+00 .COOE+00 -.650E+02 .000E+00 340E+02 .000E+00 .000E+00 .000E+00 
3 .OOOE+00 .0OOE+00 .000E+00 .0O0E+00 .000E+00 .000E100 .QOOE+00 .000E+00 .000E+00 

Element fixities 

el. fix-x fix-y fix-r x_local y_local 

1 1 1 1 .000E+00 .000E+00 

2 0 0 0 .000E+00 .000E100 

3 0 0 O .OOOE+00 .000E+00 


Description of Parameter Names in Element Data Structure 
ILD number Description 
1 SMALLEST POINTER 
2 ELEMENT SHAPE(E.G. SPHERICAL) 
3 X-~-LOCATION OF ELEMENT CENTROID 
4 Y-LOCATION OF ELEMENT CENTROID 
5 ELEMENT GLOBAL ROTATION - THETA 
6 X- VELOCITY OF ELEMENT CENTROID 
7 Y-VELOCITY OF ELEMENT CENTROID 
8 ANGULAR VELOCITY OF ELEMENT - OMEGA 
9 ELEMENT LENGTH X-DIR. 
10 ELEMENT LENGTH Y-DIR. 
11 POWER OF SUPERQUAD 
12) NUMBER OF NODAL POINTS ON SUPERQUAD 
13. ATTACHMENT LENGTH GENCE ONLY) = 0.0 FOR 
14. ELEMENT DENSITY 
15 ELEMENT MASS 
16 ELEMENT MASS MOMENT OF INERTIA 
17 EXTERNAL CENTROIDAL FORCE X-DIR. 
18 EXTERNAL CENTROIDAL FORCE Y-DIR. 
19 EXTERNAL CENTROIDAL MOMENT 
20 CENTROIDAL BOUNDARY CONSTRAINT, X-DIR. 
21 CENTROIDAL BOUNDARY CONSTRAINT, Y-DIR. 
22 CENTROIDAL BOUNDARY CONSTRAINT, ROTATION 
23 TOTAL APPLIED ELEMENT LOAD, X-DIR. 
24 TOTAL APPLIED ELEMENT LOAD, Y-DIR. 
25 TOTAL APPLIED ELEMENT LOAD, ROTATION 
26 OLD ELEMENT VELOCITY AT CENTROD X-DIR. 
27 OLD ELEMENT VELOCITY AT CENTROB Y-DIR. 
28 OLD ELEMENT VELOCITY AT CENTROID ROT. 
29 NORMAL CONTACT STIFFNESS 
30 SHEAR CONTACT STIFFNESS 
31 MASS DAMPING FOR TRANSLATION 
32 MASS DAMPING FOR ROTATION 
33 LOCAL X-COORDINATE OF A FIXED PT. ON ELE 
34. LOCAL Y-COORDINATE OF A FIXED PT. ON ELE 
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35 ELEMENT KINETIC ENERGY (TRANSLATIONAL) 


36 ELEMENT KINETIC ENERGY (ROTATIONAL) 
37 ELEMENT KINETIC ENERGY (TOTAL) 

38 ELEMENT GRAVITATIONAL ENERGY 

39 NODAL SPRING ENERGY (AXIAL) 

40 NODAL SPRING ENERGY (SHEAR) 

41 NODAL SPRING ENERGY (BENDING) 


42 INCREMENTAL SPRING DISPLACEMENT (AXIAL) 
43. INCREMENTAL SPRING DISPLACEMENT (SHEAR) 


44 INCREMENTAL SPRING ROTATION 


45 ELEMENT POTENTIAL ENERGY (GRAVIT. & SPRI 


46 TOP CABLE LENGTH 

47 BOTTOM CABLE LENGTH 

48 TENSION IN TOP CABLE 

49 X-COMPONENT OF TENSION IN TOP CABLE 
50 Y-COMPONENT OF TENSION IN TOP CABLE 
51 FORCE IN TOP JOINT OF ATTENUATOR 

52. TENSION IN BOTTOM CABLE 


53 X-COMPONENT OF TENSION IN BOTTOM CABLE 
54 ‘Y-COMPONENT OF TENSION IN BOTTOM CABLE 


55 STIFFNESS OF TOP CABLE 

56 STIFFNESS OF BOTTOM CABLE 

57 INITIAL LENGTH OF BOTTOM CABLE 

58 BREAKING STRENGTH OF TOP CABLE 

59 BREAKING STRENGTH OF BOTTOM CABLE 
60 TOP CABLE PLASTIC DEFORMATION (FT) 


61 LARGEST POINTER (CHECK WITH NPRMAX !) 


Fence Support Data Parameters 

Horizontal cable span = 55.0 
Vertical cable sag = 3.50 
Horizontal distance between columns = 2.50 
Youngs modulus of cable = .200E+08 


y-coordinate of top cable connection = 6.67 
y-coordinate of bottom cable connection = -4.00 
No. of top/bottom connections on fence = 2 


No. of tire columns in fence layer = 13 
coefficient of restitution = .400 
x-section area of top cable = 1.08 
x-section area of bottom cable = 250 
coefficient of friction = .000 

slack in between bottom connections = .500 
top cable breaking strength = .100E+09 
bottom cable breaking strength = ,100E+09 
No. of fence layers = 2 


gap between fence layers = 2.60 
Maximum overlap (rock/fence) = ,6000E+00 
Maximum overlap (fence/fence) = ,5000E+00 


Element I.Ds of attenuators in fence layer: 1 
[LDnos.= 3456789 
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Element (Ds of attenuators in fence layer: 2 
LD nos. = 11 12 13 14 15 16 17 


element mass moi. 
1 .4000E+02 .1337E+04 
2 .3723E+102 .2148E+02 
3. =,2985E+02 .1708E+03 
4 .2985E+02 .1708E+03 
5 .2985E+02 .1708E+03 
6  .2985E+02 .1708E+03 
7  .2985E+02 .1708E+03 
8 .2985E+02 .1708E+03 
9 .2985E+02 .1708E+03 
10 = .2985E+02_ .1708E+03 
11 = .2985E+02 .1708E+03 
12 = .2985E+02 .1708E+03 
13. .2985E+02 .1708E+03 
14 = .2985E+02 .1708E+03 
15 = .2985E+02 .1708E+03 
16 = .2985E+02 .1708E+03 
17s .2985E+02_ .1708E+03 
18 .2985E+02 .1708E+03 


Normal contact stiffness (rock/fence)= .1944E+06 
Maximum fence penetration = .5000E+00 


Normal contact stiffness (fence/fence)= .2522E+06 
Maximum fence penetration = .5000E+00 


Normal contact stiffness (ground/rock)= .8723E+08 
Information: S/R chekdt- time step= .827E-04 
note, max. stable time step = .827E-04 


max time step governed by axial stiffnesses 
ofelementno. = 3 


time=  .496E-02 secs 
el x y theta vx vy ang-v 


1 .000E+00 -.500E+01 .000E+00 .Q000E+00 .000E+00 .000E+00 
2 .216E+01 -.402E-03 .165E+00 -.490E+02 -.158E+00 .306E+02 
3 -.428E-01 .166E-05 .591E-03 -.200E+02 .179E-02 .378E+00 
4 -.228E-04 .545E-01 .266E-04 -.204E-01 -.215E-02 .237E-01 
5 -.429E-07 .218E100 .500E-07 -.517E-04 -.148E-05 .603E-04 
6 -.474E-10 490E+00 .552E-10 -.725E-07 -.149E-08 .845E-07 
7 -336E-13 .871E+00 .392E-13 -.627E-10 -.201E-11 .730E-10 
8 -165E-16 .136E+01 .192E-16 -.361E-13 .148E-12 .421E-13 
9 -.585E-20 .196E+01 .682E-20 -.147E-16 -.280E-12 .172E-16 
10 .0O0OE+00 .550E+01 .000E+00 .000E+00 .000E+00 .000E+00 
11 -.260E+01 .384E-17 .OOOE+00 .QOOE+00 .152E-14 .000E+00 
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12 -.260E+01 
13 -.260E+01 
14 -.260E+01 
15 -.260E+01 
16 -.260E+01 
17 -.260E+01 
18 -.260E+01 


45E-01 .000E+00 .O00E+00 .182E-15 .000E+00 
.218E+00 .000E+00 .Q00E+00 ~.109E-14 .000E+00 
490E+00 .000E+00 .000E+00 .200E-15 .000E+00 
871E+00 .000E+00 .000E+00 -.101E-14 .000E+00 
.136E+01 .000E+00 .000E+00 .148E-12 .000E+00 
-196E+01 .OOOE+00 .000E+00 -.280E-12 .000E+00 
SSOE+01 .OO0E+00 .000E+00 .000E+00 .000E+00 


time= .993E-02 secs 


el x y 


-.194E+00 
-.791E-03 . 
-.430E-05 . 
-.156E-07 . 
-.393E-10 . 
-.711E-13 . 
-.962E-16 . 
10 .000E+00 
11 -.260E+01 
12 -.260E+01 
13 -.260E+01 
14 -.260E+01 
15 -.260E+01 
16 -.260E+01 
17 -.260E+01 
18 -.260E+01 


OMO~ IDM PWYH 


theta vx vy ang-v 


.000E+00 -.500E+01 .000E+00 .000E+00 .000E+00 .000E+00 
.196E+01 -.154E-02 .310E+00 -.332E+02 -.295E+00 .285E+02 


.262E-03 .160E-02 -.377E+02 .224E+00 -.945E+00 
543E-01 .922E-03 -.454E+00 -.988E-01 .529E+00 
218E+00 .SOIE-O5 -.303E-02 -.461E-03 .354E-02 
490E+00 .182E-07 -.134E-04 -.167E-05 .156E-04 
871E+00 .458E-10 -.400E-07 -.799E-08 .466E-07 
136E+01 .829E-13 -.843E-10 -.449E-10 .982E-10 
196E+01 .112E-15 -.130E-12 -.778E-12 .152E-12 
S50E+01 .0O0E+00 .000E+00 .OO0E+00 .000E+00 
-.123E-04 .139E-03 -.110E+01 -.557E-01 .537E+00 
545E-01 .562E-08 -.353E-04 -.858E-06 .411E-04 
.218E+00 .679E-13 -.547E-09 -.140E-10 .638E-09 
.490E+00 .332E-18 -.312E-14 379E-15 .363E-14 
.871E+00 .0O0E+00 .000E+00 -.202E-14 .000E+00 
.136E+01 .000E+00 .000E+00 .296E-12 .0O0E+00 
.196E+01 .000E+00 .000E+00 -.561E-12 .000E+00 
SSOE+t01 .0O0OE+00 .000E+00 .000E+00 .000E+00 


time= .149E-01 secs 


el x y 


1 .G00E+00 - 
2 .183E+01 - 
3 -.357E+00 
4 -.725E-02 
5 -.940E-04 
6 -.638E-06 
7 -.314E-08 
8 -.116E-10 . 
9 -331E-13 
10 .000E+00 
11 -.266E+01 
12 -.260E+01 
13 ~.260E+01 
14 -.260E+01 
15 -.260E+01 
16 -.260E+01 
17 -.260E+01 
18 -.260E+01 


theta vx vy ang-v 


SO0E+01 .000E+00 .000E+00 .000E+00 .000E+00 
342E-02 444E+00 -~.211E+02 --477E+00 .252E+02 
-769E-03 -.800E-02 -.307E+02 -.770E+00 -.344E+01 


530E-01 .845E-02 -.225E+01 -.405E+00 .262E+01 
.218E+00 .LIQE-03 -.494E-01 -.248E-O01 .576E-01 
A90E+00 .744E-06 -.399E-03 -.153E-03 .465E-03 
.871E+00 .366E-08 -.229E-05 -.133E-05 .267E-05 


136E+01 .135E-10 -.971E-08 -.149E-07 .113E-07 


.I96E+01 .385E-13 -.313E-10 -.198E-09 .365E-10 


SSOEH)1 .000E+00 .000E+00 .000E+00 .000E+00 
-193E-02 -.599E-03 -.194E+02 .163E+01 -.242E+00 
45E-01 .630E-04 -.363E-01 -.414E-02 .423E-01 
218E+00 .184E-06 -.L57E-03 -.769E-05 .183E-03 
ASOE+00 .281E-09 -.321E-06 -.163E-07 .374E-06 
.871E+00 .260E-12 -.370E-09 -.397E-10 .431E-09 
.136E+01 .160E-15 -.272E-12 .357E-12 .317E-12 
-196E+01 .408E-19 -.109E-15 -.841E-12 .126E-15 
SS0E+01 .00OE+00 .000E+00 .000E+00 .000E+00 
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IL5.3 Sample Geometry Plotting Session 


The following output is from the computer screen during a sample session using GPLOT, 
the interactive geometry plotting module. This session produces an animation of the 
sample2 analysis: 


NOTES: 
(i) The command g executes a batch file g.bat which executes gplot from a subdirectory. 
(ii) The user responses are prefixed with >. 


78 
Enter input-data file name (Exclude Ext.) : sample2 
Enter Plotting Time: 
>11 
do you want to plot element numbers (y/n)? 
mail 
current window coordinates are: 
xmin y-min x-max y-max 
-100 -150 100 9.69 
do you want to change the window coords (y/n)? 
> 
animate or single frame (a/s) 7 
>a 


Stop - Program terminated. 


The hardcopy output from grock.ps is shown in Figure II.5 below. 


ROCK WEIGHT = 1282 LBS 


Figure II.5 A Geometry Snapshot Plot 
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X - VELOCITY OF ELEMENT CENTROID 


IL5.4 Sample Time History Plotting Session 


The following output is from the computer screen during a sample session using HPLOT 
the interactive time history plotting module. This session produces a time history of the x- 
component of the rock velocity throughout the impact for the sample2 analysis. A 
hardcopy file called vrock.ps is created during the session. The hardcopy output from 
vrock.ps is shown in Figure II.6: 


NOTE: The command h executes a batch file h.bat which executes hplot from a 
subdirectory. 


ROCK WEIGHT = 1202 LBS 


ig! 
a5 


TIME. IN: (SECS33 


Figure II.6 Time History of x-component of Rock Velocity 
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C: \GRAHAM\FENCE oh 


C: \GRAHAM\FENCE *hoplat\hist 
ENTER HISTORY PLOT FILENAME (26 CHARACTERS MAX.) 


sample2 

Select from Available Elements Numbers: 
2 3 4 Ss & 

Enter Element Number: 

2 


Element Number 2 is Selected 
Do you want to select another element? 
Enter yes) or nto) 
n 
Select y-values from available components: 
Case: 1 FORCE IN TOP JOINT OF ATTENUATOR 
Case: 2 TENSION IN BOTTOM CABLE 
Case: 3 TENSIGN IN TOF CABLE 
Case: 4 X — VELOCITY OF ELEMENT CENTROID 
Case: 3S Y - VELOCITY OF ELEMENT CENTROID 
Cases 6 ANGULAR VELOCITY OF ELEMENT —- OMEGA 
Case: 7 ELEMENT KINETIC ENERGY (TRANSLATIGNAL) 
Case: 8 ELEMENT KINETIC ENERGY (ROTATIONAL) 
Enter Case Number: 
4 


Case: 7 ELEMENT KINETIC ENERGY (TRANSLATIONAL) 
Case:, 8 ELEMENT KINETIC ENERGY (ROTATIONAL) 
Enter Case Number: 

4 

X ~— VELOCITY OF ELEMENT CENTROID is Selected 
Do you want to plot element x-values 
n 

Enter Flotting Time: 

il 

do you want to change plotting titles ? 

ROCK WEIGHT = 1200 LAS 


n 

do you want to change plotting x-title ? 

TIME IN (SECS.) 

n 

Go you want to change plotting y-title ? 
X -— VELOCITY OF ELEMENT CENTROID 

n 

do you want to change ymin, ymax ? 

ymin = -.6500E+02 ymax = - 7456E+01 

hn . 

do you want to change no. of x and y grid divisions ? 

ndivx = wo ndivy = 2 

n 
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fFLLECECEECCOL ELE LE CE COE CLL ELE LEC EE LECCE ECL EL CO ECEL EEC EE CEL OE LOE 
CEOCEECELO COO CE ECC LEC ECE EOC ELEC OCC EEOC ELE ECL OL LEE COE Eee 
COLOCOE CECE LL EEC LEE CEE LOLOL CECE CELL EC ECE CL EL CECE ECE LL CEE ECE OE 
CECELCLECECLE CEC COOL ECOL ELE COCO CCE CELLO CL ECOL CEO LL OOO CEL OELC ELIE 
LECEELCO CLEC CEL EEOC CELLO COL ELLE CLL EE EE CEO OL OL LO CE ELEC EOC EE CELE LE 
CLLCOCCOCELCOCCELCLOCOCOOE LCCC CEL EC COOOL CO LCL EL COCO LEE COOOL OCLC ELE 
ECLECLLOEL CEC EL OCLC EOE ECE ELE CCL ELLE OLE EEE ELLE LEC ELC CEL LEE EE 
CLCCCCLLL ELLE ASOREEN CIMAGE OFC T4MEL HISTORY Y LLL ELLE CLL ELE LCE 
ELEEEL ECE COOOL ELE CCL E EEC EC ECE ECE ELEC CCE LCE LCE CECE LCL ELE CELL LEE 
COCOECCLCECECOCCEC CCE CEC COC OCOCECELC ECC CLE CE LLL OLE CLE CEL CELE CEE EOE 
LCOCCECCEECEL LCE CLEC ELE CO OOO CL OL ECL OLE LCE COE LCE EEOC LE LEC LE LE CECE LE 
CLOEECLCECECEC OL LEC COCL OL COCO CEC LEE CEC EL COL EEL CLO CO CLC L ELL ELLE 


_ DO you want a hardcopy output file (y/nm)? 


y 
enter filename —- max 20 characters 


vrock.ps 
Hardcopy file is vrock.ps 


do you want to dump history data to a file (y/n) ? 


y 

enter data file name, 40 characters max. 
vrock 

data file name is: 

vrock. RAW 


enter return to continue 


Do you, want to do another plot,new plot file or quit? 
Enter -a (nother ) /n (ew) /q (uit) . 


Stop - Program terminated. 
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